#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <getopt.h>
#include <sys/time.h>
#include <gsl/gsl_cdf.h>
#include "tree.h"

int main(int argc, char *argv[]) {
    treePars p;
    treeResults r;
    struct timeval start, end;
    
    getPars(argc, argv, &p);
    
    gettimeofday(&start, NULL);
    binomtree(&p, &r);
    gettimeofday(&end, NULL);
    
    double d1, d2, c, delta;

    d1 = (log(p.S0 / p.strikePrice) + (p.r + p.sigma*p.sigma * 0.5) * p.endTime)
                / (p.sigma * sqrt(p.endTime));
    d2 = d1 - p.sigma * sqrt(p.endTime);
    c = p.S0 * gsl_cdf_ugaussian_P(d1) - p.strikePrice * exp(-p.r * p.endTime) * gsl_cdf_ugaussian_P(d2);
    delta = gsl_cdf_ugaussian_P(d1);
    
    printf("%lf %lf    ", r.optionValue - c, r.hedge - delta);
    printf("%.3lf seconds\n",\
            end.tv_sec - start.tv_sec +  (end.tv_usec - start.tv_usec) * 1.e-6);
    return 0;
}

void showtree2(int n, double **t) {
    int i, j;
    for(j = 0; j < n; j++) {
        for(i = 0; i < j; i++) printf("        ");
        for(i = j; i < n; i++) {
            printf("%8.2f", t[i][j]);
        }
        printf("\n");
    }
}

void showtree(int n, double **t) {
    int i, j;
    int srow, row, col;
    int inbot=0;
    for(i=0; i<2*n-1; i++) {
        if(i>=n) inbot = 1;
        srow = abs(n-i-1);
        for(j=0; j < srow; j++) {
            printf("       ");
        }
        col = inbot*srow;
        row = srow;
        for(j=0; j < (n+1-srow)/2; j++) {   
            printf("%7.2f", t[row][col]);
            printf("       ");
            row += 2;
            col++;
        }
        printf("\n");
    }
}

void binomtree(treePars *pars, treeResults *results) {
    int i, j;
    double **f, **S;            // f will store the option values, S the stock prize
    double optionValue;         // Will store the computed option value
    double hedge;               // Will store the computed hedge parameter
    
    int N = pars->nsteps;                           // The number of steps
    double deltaT = pars->endTime / N;              // The timestep
    double u = exp(pars->sigma * sqrt(deltaT));     // Increase factor per timestep
    double d = exp(-pars->sigma * sqrt(deltaT));    // Decrease factor per timestep
    double a = exp(pars->r * deltaT);               // Interest factor per timestep
    double p = (a - d) / (2. * sinh(pars->sigma * sqrt(deltaT)));   // Probability of stock going up
    double S0 = pars->S0;                           // Initial stock price
    double strikePrice = pars->strikePrice;
    p = p > 1. ? 1. : p;
    
    // Allocate space for the matrices, which are triangular
    f = (double **) malloc((N+1) * sizeof(double *));
    S = (double **) malloc((N+1) * sizeof(double *));
    
    for(i=0; i < N+1; i++) {
        f[i] = (double *) malloc((i+1) * sizeof(double));
        S[i] = (double *) malloc((i+1) * sizeof(double));
    }
    
    // Initialize the stock price and compute stock prices at later times
    S[0][0] = S0;
    for(i = 1; i < N+1; i++) {
        for(j = 0; j < i; j++) {
            S[i][j] = u * S[i-1][j];
        }
        S[i][i] = d * S[i-1][i-1];
    }

    // Calculate f from back to front, first set the last row
    for(j = 0; j < N+1; j++) {
        if(pars->isCall)
            f[N][j] = doublemax(S[N][j] - strikePrice, 0.0);    // Call option
        else
            f[N][j] = doublemax(strikePrice - S[N][j], 0.0);    // Put option
    }
    
    // Now recursively go back in the tree
    for(i = N-1; i >= 0; i--) {
        for(j=0; j < i+1; j++) {
            f[i][j] = p * f[i+1][j] + (1.0 - p) * f[i+1][j+1];
            f[i][j] /= a;
            if(!pars->isEuropean) {
                if(pars->isCall) {
                    f[i][j] = doublemax(f[i][j], S[i][j] - strikePrice);    // American call option
                } else {
                    f[i][j] = doublemax(f[i][j], strikePrice - S[i][j]);    // American put option
                }
            }
        }
    }
    optionValue = f[0][0];
    hedge = (f[1][0] - f[1][1]) / ((u - d) * S0);
    
    // Output the tree
    if(pars->doPrint) {
        printf("Stock prices in the tree\n");
        showtree(N+1, S);
        
        printf("\nOption values in the tree\n");
        showtree(N+1, f);
    }
    
    for(i=0; i < N+1; i++) {
        free(f[i]);
        free(S[i]);
    }
    free(f);
    free(S);
    
    results->optionValue = optionValue;
    results->hedge = hedge;
}

void getPars(int argc, char *argv[], treePars *pars) {
    pars->isCall = 0;
    pars->isEuropean = 0;
    pars->doPrint = 0;
    pars->endTime = 1.;
    pars->nsteps = 500;
    pars->strikePrice = 99.;
    pars->S0 = 100.;
    pars->sigma = 0.2;
    pars->r = 0.06;
    
    int ch;
    while ((ch = getopt(argc, argv, "cEvn:t:k:S:s:r:")) != -1) {
        switch(ch) {
        case 'c': pars->isCall = 1; break;
        case 'E': pars->isEuropean = 1; break;
        case 'v': pars->doPrint = 1; break;
        case 'n': pars->nsteps = strtol(optarg, 0, 10); break;
        case 't': pars->endTime = strtod(optarg, 0); break;
        case 'k': pars->strikePrice = strtod(optarg, 0); break;
        case 'S': pars->S0 = strtod(optarg, 0); break;
        case 's': pars->sigma = strtod(optarg, 0); break;
        case 'r': pars->r = strtod(optarg, 0); break;
        }
    }
    
    printf("Parameters:\
            \n isCall (-c)     \t%d\
            \n isEuropean (-E) \t%d\
            \n doPrint (-v)    \t%d\
            \n nsteps (-n)     \t%d\
            \n endTime (-t)    \t%.2lf\
            \n strikePrice (-k)\t%.2lf\
            \n S0 (-S)         \t%.2lf\
            \n sigma (-s)      \t%.2lf\
            \n r (-r)          \t%.2lf\n\n",
            pars->isCall,
            pars->isEuropean,
            pars->doPrint,
            pars->nsteps,
            pars->endTime,
            pars->strikePrice,
            pars->S0,
            pars->sigma,
            pars->r );
}

